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Abstract 

A Galerkin method is developed to solve the time-dependent Dirac equation in prolate spheroidal coordinates 
for an electron-molecular two-center system. The initial state is evaluated from a variational principle using a 
kinetic/atomic balanced basis, which allows for an efficient and accurate determination of the Dirac spectrum 
and eigenfunctions. B-spline basis functions are used to obtain high accuracy. This numerical method is used 
to compute the energy spectrum of the two-center problem and then the evolution of eigenstate wavefunctions 
in an external electromagnetic field. 

Keywords: Dirac equation, prolate spheroidal coordinates, two-center system, Galerkin method, 
variational method, B-spline basis set, atomic/kinetic balance 


1. Introduction 

In the last few decades, there has been a surge of interest for the numerical solution of the Dirac equation in 
many areas of physics and chemistry, motivated mostly by new advances in computational architectures and 
numerical methods, which allow to tackle complex physical problems. One specihc field that has benefited 
from these advances is laser-matter interaction where it is now possible to reach laser intensities of 10^° 
W/cm^ [1] and higher in laboratories. The theoretical description of matter subject to such intense radiation 
can only be described by relativistic quantum mechanics which requires solutions of the Dirac equation [2]. 
Traditionally however, the Dirac equation has been studied mostly in the context of relativistic heavy ion 
collisions, where the search for positron production from Uranium nuclear collisions is one of the main 
impetus [3]. 

Various numerical methods have been developed to solve the relativistic equation as analytical approaches 
are often challenging and only perturbative. Among the most popular approach is the operator splitting 
method, where the Dirac operator is separated into a set of simpler equations. Each of these resulting 
equations can then be solved by resorting to well-known and accurate numerical schemes. As there exists 
many possible decompositions of the Dirac operator, there also exists many variations of the operator splitting 
method. It is often combined with spectral methods whereby the kinetic operator is solved by the Fourier 
Transform methods while the mass and potential terms, being local operators in “real space”, can be dealt 
with by accurate approximations of time-ordered exponentials. This technique has been used in misiisiiz] 
for the Dirac equation and in mm for the coupled Maxwell-Dirac equation that includes the interaction and 
the backreaction on the electromagnetic field. Another possible decomposition of the Dirac Hamiltonian was 
given in uniEiiii] using Alternate Direction Iteration. In this case, the spin is kept aligned with the direction 
of propagation at each step (using a specific rotation in spinor space) such that simple analytical solutions 
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can be found using the method of characteristics. The resulting scheme, sometimes called “Quantum Lattice 
Boltzmann”, can be parallelized very efficiently. It can also be adapted to treat the cylindrical coordinate 
case [13] and nonlinear Dirac equations [14] . 

Although these approaches are very powerful and have very interesting properties, they are inefficient 
for finding the initial state of the system in a confining potential with bound and continuum states. Within 
the operator splitting method, these states are usually determined from a relativistic variant of the Feit- 
Fleck method 111311 US). The latter allows for the computation of the spectrum and the determination of 
bound states from a filtering technique on the time evolution of the wavefunction, thus having a very slow 
convergence. Therefore, these methods are impractical for problems in Quantum Electrodynamics requiring 
sums over all states of the spectrum. For these reasons, other approaches have been considered. One 
possibility is the use of the mapped Fourier grid, which allows to evaluate both the spectrum and the time 
evolution of the wavefunction [13 El- One problem however with this scheme is the appearance of spurious 
states, which are unphysical states created during the discrete evolution process. Direct approaches, where 
the Dirac operator is discretized without splitting, have also been attempted. For instance, implicit finite 
difference schemes can be found in [nnnmiinj while an explicit scheme is in |21] . 

Conversely, there exists very powerful schemes to solve the time-independent Dirac equation, based 
on variational methods and basis set expansion. The most important issue in this case is the variational 
collapse [23 , which is related to the fact that the spectrum of the Dirac equation is not bounded from below 
(or above). This induces spurious states in the spectrum obtained from the usual Rayleigh-Ritz variational 
method. This phenomenon is also called spectral pollution [23|- There has been several (successful) attempts 
to solve this problem and there now exists two main lines of development: 

1. New variational principles 

2. Balance principles 

The hrst case corresponds to a modification of the usual Rayleigh-Ritz minmax principle. This was first 
investigated by Talman [24] and was generalized in [23 126] • This has led to numerical methods free of 
spurious states, but which requires the solution of a nonlinear eigenvalue problem (see [22] for instance). 
The latter usually requires an iteration method and thus, necessitates a lot of computation time. The second 
case corresponds to a modification of the basis function expansion such that spinor components are related 
in some ways. This was first introduced as an empirical rule to get rid of spurious states |22j and was then 
analyzed by comparing with the non-relativistic results |281129] . However, the rigorous analysis of these 
methods is fairly recent |23| . There exists three well-known variations of the balance principle: 

1. Kinetically balanced basis function [5D] 

2. Atomic balanced basis function m 

3. Dual kinetic balanced basis function [32] 

In each of them, a different relation is imposed between basis functions of the large and small spinor compo¬ 
nents. In this work, the atomic as well as kinetic balance will be used to compute the initial state (Cauchy 
data) for the time-dependent Galerkin method. 

The Galerkin method has been applied to the Dirac equation using different coordinate systems and basis 
sets [231311 [S3 > both time-dependent and time-independent cases. In this article, we develop numerical 
schemes to study the two-center problem in an external electromagnetic field. This system has also been 
investigated extensively, mostly in connection with heavy ion collisions and heavy ion spectroscopy. The 
static case can be found in [13133113133133130113111321133133] ( an analytical approximation can be found 
in [33 ) but less is known for the dynamic case II3I33I2I1I33- The main goal of this article is to give a 
variant of these methods, based on atomic balance and B-spline basis sets. 

This article is separated as follows. In Section [^ we describe the Dirac equation studied in this article. 
Section is devoted to the derivation and analysis of the Galerkin solver for the time-independent Dirac 
Hamiltonian. In Section we derive from the Time Independent Dirac Equation (TIDE) solver, a Time 
Dependent Dirac Equation (TDDE) version. Some mathematical properties of the derived schemes are also 
proposed in this section. Some important details of the numerical implementation are given in Section [3 
along with some performance benchmarks. The numerical results are presented for TIDE and TDDE in 
Section]^ We finally conclude in Section]^ 
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2. Dirac Equation 

The Dirac equation is a quantum wave equation that describes the relativistic dynamics of spin-f particles 
(fermions) such as the electron. In this setting, the particle under consideration is characterized by a four- 
component spinor 

= [<P,xf & C\0,T;L^R^&)), 

for some positive time T. The bispinors (j)^x ^ (O, T; C^)) are usually, respectively called the 

large and small components. In the time-independent case, where we consider an interaction with a nucleus 
defined by a static external Coulomb potential 14, the wavefunction obeys the following Dirac equatior0 

idt'^ = i?o^, with Hq = ccx ■ Tp + me?(3 + I 4 (x)l 4 , 


where o: = [ax,Oiy,az) are the Dirac matrices, Hq is the Hamiltonian operator, p = —zV is the momentum 
operator, c is the light velocity, m is the electron mass, and 'I' is the four component spinor. The matrix 
structure is given by a and /3 in M 4 (C): 


CXi — 


0 (Ti 

a, 0 


and 


I 2 0 

0 —12 


where Ui are the usual Pauli matrices. The latter are 



■ 0 

1 ■ 


■ 0 -i ' 


'I O' 

- 

I 

0 

1 ~ 

i 0 

and az = 

0 -1 


( 1 ) 

( 2 ) 


We consider now the relativistic spin-f quantum particle subject to a classical electromagnetic field 
(A, V) G X K_|_, K'^). We will assume that the electromagnetic field is given at any time and that the 

back-reaction of the particle on the electromagnetic field is neglected. Therefore, Maxwell’s equations are 
not solved numerically and we parametrize the electromagnetic field by an analytical form, given below (we 
refer to [iiiiniii] for a full Maxwell-Dirac equation solver based on another approach). The equation we 
consider is then: 


idt^ = , H = a. ■ (y — ieV — eA) -|- me?(3 + (I4(a^) -I- V(t, a;))l 4 , 

where the electromagnetic field was added by the minimal coupling prescription, which guarantees a gauge 
invariant formulation. In explicit calculations however, a specific gauge is chosen: we choose the Coulomb 
gauge V • A = 0 such that the Coulomb law can be used to describe the static charged nuclei. We also set 
V = Q such that the laser field is characterized by the vectorial potential. 

This last equation gives a consistent description of bound electrons in molecules in the Born-Oppenheimer 
approximation, i.e. when the nuclei are fixed in space and included in the potential term 14- This is a valid 
approximation when the mass of the nucleus is much larger than the mass of the electron thus neglecting 
momentum exchange between photons, electrons and nuclei, which will always be the case for the systems 
considered in this study. 


2.1. Dirac equation for the two-center system, prolate spheroidal coordinates and boundary conditions 

We focus on the simple electron molecular two-center system where we consider two nuclei described by 
the Coulomb potential, as 


K = - 


Zie 


Z^e 


y/x'^ + 3- {z — Ry ^x'^ -\- {z -\- R? 


(3) 


^All the calculations will be performed in atomic units (a.u.) where m = 1, ?l = 1 and c = 1/a where we take a ^ 
1/137.035999679 as the fine structure constant. In all the equations however, we are keeping the mass explicitly, allowing to 
switch easily from atomic to natural units. 
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where Zi _2 are the nuclear charges, R is the internuclear distance and x,y,z are Cartesian coordinates. 

To treat this system, it is convenient first to consider cylindrical coordinates where 

X = r cos{0), y = r sm{9), (4) 

where r = x'^ + y'^ is the radial distance and 9 = tan~^{y/x) is the azimuthal angle. Assuming that the 

Dirac equation has an azimuthal symmetry, which occurs when the electrodynamic potential does not depend 
on 0, it is then possible to reduce the number of dimensions from 3 to 2 by separation of variables. The 
0-dependence can be factorized by using the following ansatz for the four-spinor with cylindrical symmetry 
[471I39]: 


4'(x,t) 


Mt,r, 


(5) 


where fj,i ^2 ■= jz T 1/2 and where jz is the angular momentum projection on the 2 ;-axis (it can take one of 
the values jz = ■ ■ ■ Substituting in the Dirac equation leads to 


idt'ip(t,r,z) = 

-\-CXz 


—icdr — ic- - eArit, r, z) 

2r 


c— - eAg{t,r,z) 


-icdz - eAz{t,r,z) 


+ j3mc^ -\- eVc{r, z) r, z). 


( 6 ) 


Then, by using the symmetry of the coordinate transformation and by assuming that the wave function is 
regular enough, it is demonstrated in m that the wave function can be written as 


z) 

= rl'"il(^i(<,r^,z). 

(7) 

'ip 2 {t,r, z) 


(8) 

z) 

= rl'"il(^3(t,r2,z). 

(9) 

ip 4 (,t,r, z) 

= rl'""l(^4(<,r2,z). 

(10) 


where (p admits a Taylor expansion in around r = 0. Therefore, the boundary conditions at r = 0 on 
'i/' is a Robin condition which depends on the value of pi and y ,2 |13j . These boundary conditions will be 
included in the numerical scheme with the addition of a prefactor in basis functions [391138] . 

For the two-center problem, it is known that prolate spheroidal coordinates yield more accurate results 
in both the relativistic and non-relativistic cases. Moreover, in these coordinates, the nuclei are positioned 
at the corners of the domain, facilitating the numerical implementation. For these reasons, we now turn to 
these coordinates. The prolate spheroidal coordinates which are related to cylindrical coordinates as follows 

r = , (11) 

z = R^V, (12) 


where ^ S [IjOo), rj G [—1,1] and 0 = [0,27r] (azimuthal angle). This choice is particularly attractive when 
dealing with a two center potential. To obtain the Dirac equation in these coordinates, one simply uses the 
mapping in Eqs. 0 and ([T^ along with the derivatives 


dr 

dz 


v/(g2-i)(i-^2) 
R (^2 _ ^ 2 ) 


- r]dn ], 


- 1) 

R{e - 0^) 


yd^ + 


R (^2 _ ^ 2 ) 


^drf. 


(13) 

(14) 
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2.2. Dirac equation in the time-independent case: Cauchy data 

The goal of this paper is to accurately solve the TDDE for particles subject to a classical electromagnetic 
field. Prior to this, we first determine the initial data of the Cauchy problem which is naturally chosen as 
the ground (or any bound) state of the Dirac Hamiltonian. We are then require to solve the TIDE: 


Ho'tjj{x) = Eijj{x). 


(15) 


It is convenient to write the four-spinor as i^ix) = [^(x), x(a:)]'^ S L^(IR^,C^) where (j){x) and x(x) are the 
large and small components, respectively. The eigenvalue problem (15) reduces explicitly to 


Vc{x) -\- m<? 

Rq 

4>{x) 

= E 

(j){x) 

Rq 

Vc(x) — mc^ 

_ X{x) , 

. X{x) _ 


(16) 


where 


Rn . — (Tt. 


—icdr — ic 


2r 


Equation (16) is equivalent to 


Rox{x) 

Ro4>{x) 


[E — me? — Vc{x)\(j){x) 
[E m(? — Vc{x)]x{x) 


(17) 


(18) 

(19) 


which is the common starting point for the numerical method that follows. The small component can then 
be written in terms of the large component yielding 


X{x) = 


Rq 


:Hx) 


E mc^ — Vc{x) 

This relation will be important for the analysis that follows concerning balance principles. 


( 20 ) 


3. Time Independent Dirac Equation Solver 

In this section, the numerical method used to compute the TIDE is described. As stated above, this is 
required to obtain the initial state of the time evolution of the wavefunction. The latter will be given in the 
next section. 

3.1. Rayleigh-Ritz method 

The Rayleigh-Ritz method is based on a variational principle which allows to estimate the eigenvalues of 
a given operator. These eigenvalues can be characterized by the following variational principle: 

_ (V'|-ffoll/’)L2(R3_C4) . 

(V'IV')l=(K3,C4) ’ ^ ^ 

which is nothing but the usual Rayleigh-Ritz coefficient. Finding the eigenvalue by this minimization pro¬ 
cedure is equivalent to finding the stationary point of the functional 

£\^\ = (^/>|77o|V^)l2(r 3^C4) — if(l/'|V')L2(R3_C4) (22) 

where the energy becomes a Lagrangian multiplier. This form will be used in the following to convert the 
basis set expansion into a generalized eigenvalue problem. It is well-known that the convergence of this 
method depends on the fact that the spectrum is bounded from below. This is not the case for the Dirac 
operator, owing to the presence of the negative energy states and this may induce spurious states in the 
spectrum. This is discussed in the next section. 
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3.2. About spectral pollution 

This section is an non-exhaustive summary of some key results about spectral pollution for the approxi¬ 
mate Dirac Hamiltonian, constructed using Galerkin’s techniques and balance principles. In this approach, 
one still uses the Rayleigh-Ritz variational principle, but with a different set of basis functions which ap¬ 


proximates the relation between small and large spinor components given in Eq. (20). 

A spurious state can be defined rigorously as follows. Notations, proofs and additional results can be 
found in |23| . We just summarize some key ideas of this very strong and quite technical work. We consider 
an operator A of domain D{A) C H, where iJ is a Hilbert space. An eigenvalue A G K is said spurious for 
Operator A if there exists a sequence of finite dimensional vector spaces (I4,)n O D{A) and Vn Q Vn+i such 
that 

• = D{A) 

• lim„ dist(A, (t(A|v-^)) =0 


A ^ o'(A) 


The last item emphasizes that the eigenvalue is not in the spectrum of the operator and thus, non-physical. 
Therefore, a strategy has to be developed to eliminate these states. One possibility is the use of balanced 
basis functions which are defined in the following way and for which we summarize some properties from 

[23]. 

The 4-component spinor is split via a projector operator P : H ^ H, defined by P[(j), x]^ = [^i 0]^ for (j) 
and X A balanced operator L : D{L) C PH —(1 — P)H is defined from P as follows 


• L is 1-1 


• D{L) 0 LD{L) is a core of A 

The notion of a spurious eigenvalue of Operator A associated to a projector P and balanced operator L can 
finally be defined. Assuming that there exists a sequence of finite dimensional vectors spaces {V.^)n such 
that V.^ C D{L) and C and 

- 7—a - 

•U„(14+0LR+) =D{A) 

• lim„dist(A,cr(A|y+g,j;^^+)) =0 

• A ^ a{A) 

The corresponding spurious spectrum is denoted Spu(A, P, L). In this framework, the kinetically balanced 
operator is defined by 


1 

Lkb = X-ya • P 

Zmc^ 

and the atomic balanced operator is defined 


Lab = 


1 


2mc^ — Vr 


at ■ p 


Two of the main theorems of |23| state 

Theorem 3.1. Assuming that 14 is of the form Vc(x) = —k\x\~^ for k G (0,3/2) (which includes Coulomb 
potentials) then 


Spu{Ho + V„P,Lkb) = [-1,1] 


6 








Theorem 3.2. Assuming that is such that Vc{x) ^ —k|x| ^ for k G (0,3/2) with sup(V'c) <2, (2 — 
Vc)~^Wc G and max(14, 0) G with p > 3 and Vc{x) -G-oc 0, then 


Spu{Ho + Vc, P, Lab^ = [-1, -1 + sup(yc)] 

In particular for Coulomb potentials, the spurious spectrum is always empty. 


According to these results for Coulomb potentials, spectral pollution can be generated with kinetically 
balanced bases, but not with atomic balanced bases. However, it should be noted that for a given basis set, 
spurious may as well not appear. As we are interested in the two-center system with Coulomb potential, from 
the spuriousity perspective it is preferably to use the atomic balance basis set. Notice that the numerical 
tests performed below have not exhibited any spurious state with the kinetically balanced operators (which 
is however not in contradiction with Theorem 3.1). 


3.3. Variational method and balanced basis set 

The basis which is chosen to expand V’ is a B-spline basis constructed as follows. First, following [38], we 
expand the small component ft as: 






n—1 


(23) 


('lo't ^12") 

where a/ ’ ' are the coefficients of the basis expansion and Bn ' p) are the basis functions, for components 
1 and 2 respectively, expressed in the prolate spheroidal coordinate system described in Eqs. (11) and 


(12). The basis function can then be written as the tensor product of B-spline functions b^{x) of order k as 




(24) 


where n = [i,j] G 1?, i G [1, and j G [1, n^]. Some properties of B-splines are recalled in the next section. 
An overall factor is used to account for angular momentum dependence (33 EHl EH] . It is defined by 




(25) 


consistent with the boundary conditions in Eqs. 0 to 

Using the atomic balance approach, the lower spinor components are then expanded as follows, in the 
atomic balance case: 


Rn 


X = 


2mc2-U, En=A^^B^n^ 
and as follows in the kinetic balance case 

Rq 


^ 2mc2 I 


xiiCv) = 

N 

IC ^ 

2mc^ — Vn 

^ n=l 

X2{Cp) = 

N 

IC ^ 

2mc^ — Vc ^ 

n=l 



f-a. - 

{ 

r . 

f ji) 

-dr+ — 

1 

L r J 


5(2) _ j. 


(26) 


(27) 


In the following, the presentation is done in the atomic balance framework. We refer to f38f or Remark 3.1 
for the kinetic balance framework. In prolate spheroidal coordinates, (26) becomes 


(28) 

( 29 ) 
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These formulae should be understood as r := r{^,r]) and z := z{^,r]) where the relations are given in Eqs. 
(111 and (12) for coordinates and in Eqs. and ( [T4| ) for derivatives. 

The very first step to solve our Cauchy problem is to determine the initial condition. In physical situations, 
it is often chosen as the ground state for the considered system of particles. We then have to solve an 
eigenvalue problem: iJoV'o = ^'oV’Oj where Hq is the field-free Dirac Hamiltonian. The variational formulation 
corresponds to finding stationary points of the functional 


^ bP ] = + ™C^)'5^>)l2(k3^c2) + (-Ro^^'Ix)l2(r3_c2) 

+ (xI^0</>)l 2(B3_C2) -I- (xK^C — WC^)x)i2(R3_C2) 
— E [(</i>|</>)L2(B3,C2) — (x|x)l2(H3^C2)] , 


(30) 


which is just an explicit way of writing the well-known Rayleigh-Ritz functional equation in Eq. 
Integration by part was used to write the second term in a convenient form. The notation (j•)L2(R3 c 2 ) 
stands for the Hermitian inner product on In the following, we define 2 operators C and S by 


C[ip] = [ [mc^+ Vc]\4>\^+ {Ro(t)\x) + {x\Ro(t>) + [yc-‘mc^]\x? 

sm = [ i0p + ixp. 

dB3 


(31) 

(32) 


Here, the product (-I-) is just the spinor producl|^ Using the atomically balanced bases, as described above, 
and finding the stationary points of S by setting 


d£ 


= 0 , 


dS 


= 0 , 


d£ 


= 0 , 


d£ 


= 0 


’ aaf 

for i S {1, • • • , N}^ we obtain the following discrete generalized eigenvalue problem: 

Ca = ESa 


(33) 


(34) 


(1) (2) 
AT ’ ’ 

(2) 

• • • , a)v 



(2) 

c) ,• • • , 

and 





r'(i) 

0 

^(3) 

l_.ll 

1^(3) 

^-"12 


r 

0 

0 

0 


0 

r'laiT 

'-^11 

r.(i) 

'-^22 

p(3)T 

'-^21 

^(3) 

'^21 

/^(2) 

'^11 

1^(3) 

^^22 

1^(2) 

,S = 

0 

0 

c(l) 

■^22 

0 

0 

c(2) 

^11 

0 

gC2) 

^12 

(35) 


r.(3)T 

'-^22 

r-(2)T 

l_.ll 

/^(2) 
^^22 J 


0 

0 

qC2)T 

^12 

1 



C = 


The elements of these matrices are defined by: 




■"22 






(36) 

(37) 


^For H a two-component spinor, it is defined as (H|H) = -h E2H2. 














and 


-^-^{drB^^)Br Y^;-r\]i 

r ^ } {2mc^ — Vcr 

[C^2^],, = Jd^xi^{dM^^)id.Bf) + idrBP){drBf) 
+ ^Bf^Bf^ + ^Bf\drBf^) 


J 

,2\„2 


L’-12 J*j 


.(3)1 


.(3)] 
■"22 . 


p(2)^p(2)l (K-mc2)c' 

+ -{drB^ )B^ /(W^T^ 

J d^x!^{d,Bl^'>)idrBf) + ^B^^\d,Bf) 

(T4 — mc^)c^ 

{2mc^ — 14)2 

_^/n p(l)lp(l)l 

r ^ ^ ^ J 2mc2 - K 

J d^x!^{dM^^)id.Bf) + {drBf^drBf) 

+ ^Bf^Bf + ^Bf\drBf) 

+ ^{drBf^)Bf^\- 


.(3)1 


^ j 2mc^ — Vc 

y'd’i|(a.B'‘’)(a,Bj^>) + ‘^B\'\a,Bf) 


-(a,B<'>)(a.Bf) + 4(4B.'‘')Br)2s4^ 


j(2) 




( 38 ) 


(39) 


(40) 


(41) 


(42) 


(43) 



(44) 

(45) 
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^ J (2mc2 - Vc)2 

J d^x!^{d,B^^'>)id,Bf) + 

+ ^ + ^Bf\drBf^) 

+ ^{drBf^)Bf 


^ j (2mc2 - v;)2 

J d^x!^id,Bl^'>)idrBf) + ^B^^\d,Bf) 


(2mc2 — T4)2 


( 46 ) 


(47) 


(48) 


These last expressions can then be rewritten in prolate spheroidal coordinates. In practice, the eigenvalue 
problem in Eq. (341 is solved by a standard eigensolver for sparse matrices. The integration measure is given 
in prolate spheroidal coordinates by 


d^x = — r]^)d^dr]d9. 


(49) 


Remark 3.1. Note that Matrices C, D and S are very similar to those obtained with kinetically balanced 
bases 1S8\/ . The difference comes from the presence in atomic balance of the Vc-term in the denominators 

-- - - the latter is absent when using a kinetically balanced basis. 

(2mc2 — VcY 

Although atomic balance is, from the variational collapse viewpoint, more attractive than kinetic balance, 
the presence of 14 (atomic balance) is seen to be source of numerical discrepancy of the overall convergence 
rate, and special treatment is then necessary to tackle this additional difficulty. 

3 . 4 . Some basic facts about B-splines 

We recall some basic facts about B-splines. We refer for instance to [IS] for details. First B-splines are 
fully determined by their order k^^ri and knot vector using the following iterative formula 

&?(^) = , bt\x) + b^Ylix) (50) 

with initial conditions 

bl{x) = 1 ioT ti ^ X < ti+i and b] = 0 otherwise (51) 

where tfs are knots coordinates. The number of knots, also referred as breaking points, at a given coordinates 
essentially determines the regularity conditions at that point: the number of knots points should be maximal 
at singular points (at the Coulomb singularity position for instance) to allow for a discontinuous-like behavior. 
As in [3S], throughout this work, the knot vectors are given by the sequences 

1 — Cl ~ “ Cfcf ^ Cfc{+1 < ■ ■ ■ < Crif-l-l = • • • = = Cmax (52) 

-1 = 771 = • • • = ?7fe^ < ?7fc^+i < • • • < 77„^+i = • • • = = 1 ( 53 ) 
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Here, n^^ri are the number of spline functions in ^ and 77 coordinates respectively. 


Considering / € C^{[a,h]) (r € N) the distance between / and the space of B-splines S'^ of degree k, 
with r < k < n, is given by 

dist(/,S'^) = inf llZ-gll < fc’’ max |a;i - a:i+ir||/(’')||oo (54) 

ges^ —fc^j$;n+l 

As a consequence, assuming that the solution to the Dirac equation is regular enough, we can expect a very 
good accuracy with high order B-splines. The methods developed here are Galerkin’s methods and therefore, 
require the numerical evaluation of several integrals. In this work, Gauss’ quadrature methods will be used 
to approximate the integrals constituting the stiffness, mass matrices. We recall that 



n 

x)dx - y^ ujif{xi) 

i=l 




22"(n!)‘^ 


2 n!( 2 n-t \z\=r 


max \f{z)\ 


with {xi)i roots of Legendre’s polynomials and its weights.. 


4. Time Dependent Dirac Equation Solver 

4.I. Finite element method for TDDE 

The Cauchy problem we now consider is: 

idtf^ = Fl-tp, {t,x) G M+ X ip{0,x) = 'tpo{x),x G (55) 


The initial data fjo is taken as a state for the time-independent Dirac operator, that is an eigenfunction 
associated with one of the eigenvalues of Hq. Now in the field dependent case the TDDE can be rewritten: 


(j){t,x) 


Vc{x) + mc^ 

R 

(t>{t,x) 

X{t,x) _ 


R 

Vc(x) — mc^ 

_ X{t,x) _ 


with R := Rq — ecr ■ A, that is 


(t>{t,x) 


Vc{x) + me? 

Rq 

(j){t,x) 

1 

0 

—ecr • A 

(t>{t,x) 

X{t,x) _ 


Rq 

Vc{x) — mc^ 

_ X{t,x) _ 


—ecr • A 

0 

_ X{t,x) _ 


Where A is the vector potential corresponding to to some field E, E = —dtA and B = V x A, in the 
Coulomb gauge, where V = Q. 

To get a Galerkin method from the preceding equation, we have to project on basis functions. To perform 
this procedure, we introduce a basis spline B defined by the atomically balanced procedure described above, 
which give the j’th basis function spinor as 


Br.= 


<(’2j 

Xi.i 

X2,i. 


ic 

2mc^ — Vc 
ic 

2 mc^ — Vc 


[-dr - f] 


(56) 


Then, the weak form of the Dirac equation is obtained as 


{Bj\idtlf)L2(M3^Cr) = J € {I,”' ,Af} 


(57) 
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where as usual, the test functions were chosen as the basis function spinor B. The last equation can be 
rewritten more explicitly as 

|ll9t(/))i2(R3_C2) + (Xil*'9tX)L2(R3,C2) = {4>j\{^c + WC^)(/))i2(R3_C2) + (Xj | (K ” )x) L2 (R3 _C2 ) 

+ ((/)j|i?0X)L2(R3_C2) + (Xi|-R0^)L2(R3_C2) 


-e{(j)j\{(T ■ A)x)i2(R3_c2) - e(xj|(cr • A)0)i2(R3_c2), (58) 

for j e {1, • • • , A}. These equations are then discretized by using a basis set expansion with time dependent 
coefficients as 

N 

n—1 
N 

(60) 

n—1 

N . . 

= 2mc^- V X!) (61) 

° n=l 1 1 

N . s 

X2it,^,v) = 2mc^- V X!) + ~ + (62) 

° n=l 1 1 

which are the atomically balanced basis functions described in the last section, but with time-dependent co¬ 
efficients. We then arrive at the semi-discrete TDDE scheme, whereby the spatial discretisation is performed, 
which writes: 

iSa(t) = (C-f D(t))a(t) 

with a(t) = [a^^^(t),--- , , ■ ■ ■ , ci^\t), c[^\t), ■ ■ ■ is the time depen¬ 

dent unknown. 

Possible time discretizations include: 


• Explicit Euler scheme, which is nonunitary: 

Sa"+i = Sa” - iAt„(C + D”)a”, 

where a" = a(t„) for n € N. 

• Semi-implicit scheme (Crank-Nicolson scheme) which is unitary: 

Sa"+i = Sa" - z^(C + D")a" - z^(C + D 


n+lNa^+l 


Sa"+i = Sa" - z^(C + D")a" - z^(C + D 


"')a"+i 


or more generally Runge-Kutta type schemes. 
Simplectic integration schemes, such as: 


i(tf) = T exp -i J S ^(C-)-D(t)) a.{U) 

= exp [—zS~^ (C -f T>{ti + St/2))] a{ti) + 0{St^) 


(63) 

(64) 
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(65) 

( 66 ) 



Matrices S, C, are identical to the ones defined in the time-independent case, in Eqs. (36l to (48|. 
The only time-dependent matrix is the one that includes the electromagnetic field D(t). It is obtained by 
discretizing in space the following terms of the weak functional: 


D[t) := -e((^j|(cr • A)x)i2(R3_c2) - e(xj|(<T • A)(/))i2(R3_c2) 
By using the basis expansion, it can be written as 


D = 


■^11 ^12 
jj(4) ]^(4) 

LU2I -'-^22 


tj(3) 

■^11 ^12 
j-)(3) j^(3) 

-'-'21 ^22 


(67) 


( 68 ) 


The entries of this matrix are 


D 


(3)] 

11 J 




D, 


(3)] 
22 J 


D 


(3)1 


D, 


(3)1 




ie J - {Ar - iAe) [-Sf ^(9,5^) + 

ie J d^xl^A,B^^\d,Bf) - (Ar + iAe) ^-B^^\drBf) - ^^B^^Bf 

ie j d^xiy-{Ar - iAe)B\^\d,Bf^) - A, [-sf ^(S.sf ^) - 
ie ( d?x\{Ar+iAe)Bf\d,Bf'^)+A, \-B^^\drBf'^) + 


2mc^ — Vc 
c 

2 mc 2 - 14 
c 

2mc^ — Vc 
c 

2 mc^ — 14, 


and 




D 


(4)1 
22 J 


D 


(4)1 

12 J 


D, 


(4)1 

21 J 


ie J d3x|-A,(a,sW)BW + (Ar+iAe) + ^B^'^Bf'^ 

ie J d^xi^-A,{d,B^^'>)Bf + (Ar - iAg) [-(^^sf ))Bf) - ^B^^^B^ 

ie j d^xiy-{Ar - iAg){d,B\^^)Bf - A, [-(5,5^)^]"^ + 

ie j d^xiy{Ar + iAg){d,Bf'>)Bf^ +A, \-{drBf'^)Bf'> - ^B^B^ 


2'mc^ — 14 

c 

2mc^ — Vc 
c 

2mc^ — Vc 
c 

2mc^ — 14 , 


(69) 

(70) 

(71) 

(72) 

(73) 

(74) 

(75) 

(76) 


Again prolate spheroidal coordinates are used to numerically evaluate these integrals. 

4-2. Mathematical properties 

The Galerkin method presented in Section |4.1| has several nice and attractive mathematical features 
which are detailed in this section. These properties are valid, except when the opposite is specified, with 
kinetically and atomically balance bases. Recall first, that the use of prolate spheroidal coordinates leads to 
a very convenient position (for local mesh refinement) of the molecule nuclei at the corners of the domain. 
The first important result is related to the structure of Matrix S (34) involved in the TDDE solver described 
in Section (4.1). 


Proposition 4.1. Matrix S defined in (34), (44) and (46) is a Hermitian matrix (real eigenvalues). 

Although in principle, 0 can be an eigenvalue, it will be necessarily unique with atomically balanced basis 
(at least), as no spectral pollution is expected in that case, |2S]. 

In the sequel, we are interested in the consistency and stability of the time dependent solver. The field-free 
TDDE: 


idt'f = Hotlj, ijj{-,0) = fo^-) 
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where i?o0o = -E’o</> 0 ) with Eq the ground state energy, has the following exact solution (/)o(-) exp(— 

When this property is satisfied at the discrete level, up to the order of the time discretization, we will say 
that the TDDE solver is consistent with the eigenvalue solver. We have 


Proposition 4.2. Assume that the time operator dt and variable, are discretized using i) an explicit Euler 
scheme or ii) a Crank-Nicolson scheme (64), then the TDDE solver (4.1) is consistent with the eigenvalue 
solver (22). 


This simple property is very important from a practical point of view. In fact, this result can be extended 
to a large class of semi-discretization in time, but for the sake of simplicity we restrict the analysis to these 
two cases. 

Proof. The numerical ground state is constructed using the same atomically balanced basis and same mesh 
as the TDDE solver. Indeed in that case, D is identically zero, and the semi-discrete scheme becomes, for 
t ^ 0 


iSa(t) = Ca(t) 

with, by assumption a(0), defined by Ca(0) = i?oSa(0). 

Using the explicit Euler scheme, we get: 

Sa^ = Sa(0) - fAtoCa(O) = Sa(0) - fAtoEoSa(O), 
which can be easily re-written as 


a^ = (l — iAtoEo)a.{0). 

By induction and for time steps At;, with I ^ 0, one obtains furthermore: 

n n 

*Sa(^At,) =Ca(^Ati) 

/=0 1=0 


Assuming that the solution at the n’th timestep is a" = njUQ^(l — iAt;Eo)a(0) and from 

Sa"+^ = (S-fAt„C)a’", 


we can obtain the solution at timestep n -I- 1 by induction: 

a”+i = (1 - iAtnEo)a.^ = niLo(l - iAtiEo)aiO), 


where a" = • • • 

result that the discretized time evolution operator is 


(l),n (2),n 


T2).n 

-'■AT 5 Cl , 


5 CjY ,Ci , 


,C 


(2),™l 

N 1- 


This leads to the expected 


n[Lo (1 - i^tiEo) = 1 - iEo ^ Ati -b 0{nAtl^) = exp{-iEo ^ Ab) -b C>(nAt^) 

1=0 1=0 

where At^o = maxo^j^n 


In the case of a (semi-implicit) Crank-Nicolson scheme, the same reasoning can be performed. For one 
time iteration, we get 


, Atn Atn , 

Sa^ = Sa(0) - i^Ca(O) - ■ 
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This is written as 


Sa^ 


Atg Atg . 

Sa(0) - i^T;gSa(0) - *^Ca\ 


from which we obtain 


(S + i^C)ai = S(1 - i^Eo)aiO). 

Then, we deduce an explicit form for the time iteration given by 

ai = (l-*^T;o)(S + *^C) 'Sa(O). 

Now for Atg small enough, this can be simplified further because 

(s + = (i - - ^(s-ic)")s-i + 0{A4). 


Then, the time evolution operator for one time iteration has the simple form 

ai ^ e"*-^*“^“a(0) + O(At^). 

From this, the discretized time evolution operator is obtained again by induction. We finally get 


= n; 


=1 


(l + f^S"^C) a(0) = e"*'®«^"=o^*'a(0)+C>(nAt3), 


z^S-^C) 

and conclude again using similar arguments as for the Euler explicit scheme. □ 


We next state some result regarding the stability of the finite element method (4.1). 


Proposition 4.3. The semi-discrete TDDE solver (4.1) with explicit Euler-based time discretization is l^- 
unstable. The semi-discrete TDDE solver with Crank-Nicolson-based time discretization, (64), is -stable. 


Proof. From the proof of Proposition |4.2[ stability is ensured in the explicit case when the spectral radius 
of the discrete evolution operator satisfies 


p(n[Lo(i-fAt,s-i(c + D'))) ^1. 


In the field-free case, the spectrum was computed using the atomically balanced method. In that case, and 
as proven in [23], there is no spurious eigenvalue and all the eigenvalues are also real. We can conclude that, 
assuming that the eigenvalue solver is exact, the explicit Euler scheme is theoretically unstable. 

In the Crank-Nicolson case, the scheme reads 


Sa”+i = Sa" - + D”)a’" - *^(C -f D"+i)a"+i 

so that, we formally have 

a"+i =n]Lg(l + i^S-i(C-|-D'+i)) D'))a° 

The requirement for stability is then that 


p(n[Lo(l + *^S-i(C-f D'+i)) '(l-z^S-i(C-f D'))) ^ I 
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We note that in the field-free case 


'(nr=o(i + *^s-ic) '(i-*^s-ic)) ^niLop((i + *^ 


S-^C) \l-i^S-^C) ) ^ 1. 


Now as S has real eigenvalues, this condition is trivially satisfied, and then as 


p((I + *^S-iC) '(l-z^S-^C)) =1 
we have, independently on the B-spline order, |a "+^|2 ^ |a°| 2 j where |a °|2 denote the ^^-norm of a°. 


At, 


In the laser-field case, we note that S^^(C -|- D") does not necessarily have real eigenvalues. By regu¬ 
larity of the electromagnetic field, we can however deduce that D"+^ = D” -f- 0(At„). We can reformulate 
the problem into 


^ I -h lAtA” ^ ^ 

^\l + iAtA^+ 0{Atl)) ^ ^ 

for some complex matrix A". The stability condition is only ensured up to a At^ term at each time iteration. 
For the same reasons as described above, the following scheme is then stable: 

Sa"+i = Sa" - ^^iC + D”)a” - i^{C + D")a"+i 


□ 

We now state an important result about the convergence of (4.1) with Crank-Nicolson semi-discrete scheme 
in time. Although a full mathematical study of the well-posedness of 


idtip = V'(0, •) ='0o(-) 


(77) 


would be necessary in order to determine the function space, the solution to (0 is living in, we can still 
give some relevant information about the convergence, without an explicit knowledge of these spaces. 


Proposition 4.4. Assume that for G H, the solution to (55), tfj, formally belongs to (0,T;V^, where 
V C L^(IR^,C^) is an Hilbert space compactly imbedded and dense in H and approximated by a finite dimen¬ 
sional vector space Vn- We also assume that {Bj)j := i<j< 7 v ® basis ofY^, such that 

—y — ^ ^ 

Vn = V. Then (4.1) with Crank-Nicolson-based time discretization is convergent. 


Sketch of the Proof. We follow the usual procedure, such as the one presented in [49] and adapting the 
proof to the Dirac case. Under the above assumptions, we define the canonical projector, Ph^, from V to 
Vn as follows 


N 

PhM'4’{tn) = ® Bj 

i=i 


with Bj € Vn and 


4’j(tn) = {'f’ifm ■)\Bj)l^(r3^c*) 


The numerical approximation is defined as follows 

N 

i=i 
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where 




and the numerical error: 


Ph^'^itn) = Ejli - V'i(in)) O Bj 


We also set: 


Pj '■— (eft„|^i)L2(B3,C'‘) 


Now from the scheme 


+ \{H{tr,+,)rht^ - i?(t«)V';:„i^,)L2(E3,c4) = o 


we get 




-^/iivV'(in)|^j)L2(K3,c4) 


^ riiv r n^aiat ^ Ai-iv r i'-'j/ 

+ - A?(t„)V'J^^|^i)L2(R3,C4) = -^{PhNi’i.tn) - ^’/ijvV'(An+l)|^j)L2(K3_C4) 

which can also be rewritten 

~ PhN'^{^n+l)\Bj) L‘^(M3 - -^{i’hM ~ PhN'‘Pi'tn)\Bj) l2(^r 3 ^c'^) 

+ 2(^(tn+l){'lPl^^ - fftjvV’(A«+l))|^i)L2(R3,C4) + 2{H{tn){'lpl^ - -P?ijv^(An)) l^i)L2(R3,C4) 

{^h]\f'4^(j^n) -^hiv V^(^'AT'+1) l■^i)I/2(E3,C^) “I” 

and becomes 

“ e^„l^i)i^(KhC4) + -^{Hitn+i)e'^'^'^\Bj)L2(R3^c*) + l^i)i='(RhC4) 

“ {^hN'4^{^n) ^^ivV’(^n+l)|^j)L2(]R3,C^) -^(^n)-P/ijv'^(^n) )l2(E^ ,C‘^) 


We set 


■ At 


-^{PhN'^i.tn) - P/i„f/'(An+l)) - l^{H{tn+l)PhM{tn+l) + H{t„)Phr,tp{tn)^ 


which is also equal to 


{£hJPj)LHR^C‘^) = {^[Phj~,^{tn) - Phj~,Htn+l)'^ 

--i7(t„+l)(P/j„'!/'(An+l)-V'(An+l,-)) - 2^itn)(^Ph«'ip{tn) - i’itn,-) 
+ ^H{tn+i)^{tn+l,-) + ^H{tn)'lp{tn,-) ^j)^2(r3 c4) 


d , , i9'i/' 

=Ph^^ 


17 



and for all j and all n ^ 1 


/ d'4’ \ 

^f/^2(R3 C4) ~ ■)l^j)L2(R3^C4) 


thus 


/ 1 / \ 1 dib 1 dib 


2 dt 

-]^H{tn+l){PhM'^{tn+l) - ^{tn+l,-)) “ (tn) V'{^n) “ V'(^r, 




L2(K3,C4) 


We now set 


and 


cn 

dh,. := 


At 


/ \ 1 / dib dib 

^\J^hM^{^n) ~ PhM'4’{^n+l)j + 2 


l/A, = 


n+1) V^(tn+1 )'0(tn+l 5 ‘)^ i,^n^ 

with Following [5^ and assuming that ijj G C^{0,T; H) we get 

At„ /■‘"+1 , 53 ^ 1 rtr.+i 

Then, we have \l3i)L^(R3^ci-^ that goes to zero when h —>■ 0, due to the density of Vn is V. 

Now, from 

- e5J„|Sj)L2(R3_c4) + -(tt(t„+i)e5^+^ + |Sj)l2(r 3^C4) = At„(d^^ + |Sj)l2(r3_c4) 

we have, without approximation 


Se 


,n+l 


Se” - i^{C + D”)e” - i^(C + D”+i)e” + At„(5” + v^) 


At 


where 


e” = - ’",••• , {tn) - a^N ^ 

• • • , • • • ,^hl(t„) - C^^’”] 

and d” = ((< 5 ;j„|^j)L 2 (R 3 ,c 4 ))^-, = ((i^^^|Sj)l 2 (r 3 _c 4 ))^.. Now we deduce 


Then 


e»+i = 


(S + *^(C + D”+i))e”+i = (S - i^(C + 
Ati, 


- At„(d” 


(d" +1/”) 


(s + *^(C + D'+i)) '(s-z^(C + D'))e" + At„(s + i^(C + D"+i)) '( 

Finally from 

\'^hM ~ ■)\h ^ - Phf,'ilj{tn)\H + |(/- P/i„)V^(t„, Olff = leljn + \{I - PhA'>Pitn,-)\H 

we formally conclude of the convergence of the method, as in [33] . We note again that this conclusion is only 
valid under strong reasonable assumptions. □ 
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5. Numerical implementation 


The numerical method described in previous sections have been implemented in a high performance 
parallel code. This is required because the calculation of physical observables entails a large amount of 
computational resources: the typical time step should obey St < 1/mc^ to guarantee high precision [12] 
while the dynamics of typical external laser fields occurs on much larger time scales. Moreover, for QED 
calculations, every negative energy states has to be evolved in time and thus, demand a large number of 
time evolution calculations. 

The parallelization is performed by using the capabilities of the PETSc [SD] and SLEPc |5T| linear algebra 
libraries. Because the B-spline basis functions have compact but overlapping support, it is not convenient to 
employ a standard domain decomposition, as in |12j , for example. Rather, the parallelization is accomplished 
by distributing the solution vector a on many processors as 


^(Olproc 1 
^(Olproc 2 


(t), af ^ (t), cf ^ (t), cf ^ (t), • • • , (t), (t ), (t), 


[«i+i W, cj^+i W, c^+i W, 


d2) 


( 1 ) 


.(2) 


( 2 )/ 


4« 


,( 2 ) 






"271 


a{t)\ 


proc M 


^N-n+l 




n+1 




iV-n+1 




AT-n+l 


it),- 


,(i), 


{t),a^^\t),c)^’(t),c%>(t)]. (78) 


Ji) 


( 2 )/ 


Here, N is the number of basis functions, M is the number of processors and n is the number of basis 
function stored on one processor: they are related hy N = M ■ n. This ordering of coefficients insures that 
all the spinor component contributions with the same support are stored on the same processor. Moreover, 
it is consistent with the PETSc parallel matrix storage, which adopts a row-wise storage type. Then, the 
number of entries for the matrices S, C and D is the same on all processors, insuring an equal load on every 
processor. The calculation of these matrices does not require any inter-processor communications and thus, 
this part of the calculation is embarrassingly parallel. Communications are required when the linear system 
or the general eigenvalue problems are solved: these operations are dealt with efficiently by the PETSc and 
SLEPc libraries. These features make for a very efficient parallel code with an excellent parallel speedup. 

In order to show the efficiency of the proposed parallelization, we study the time evolution of the wavefunction 
for a dithorium two-center system subjects to an external electric field. Data are respectively as follows: 
A(c = Njj = 10, = Njj = 20, then iVj = TV,, = 30, with B-spline order fixed to 5. The time step is fixed to 

10“® and 10^ time iterations are performed. We report in logscale in Fig. the computational time using 
respectively 1, 4, 16, 32 and 64 processors. Notice, that the discrepancy in the scalability graph, observed 
for 64 processors, is a simple consequence of the moderate size of this benchmark. 


6. Numerical results 

This section is devoted to the numerical validation of the B-spline method presented in Sections [^ 
Detailed physical properties of system under consideration, will be studied in a forthcoming paper, specifically 
dedicated to quantum relativistic particles subject to external classical field. In this paper, the considered 
particles are dihydrogen (.^ 1,2 = 1) or dithorium (.^ 1,2 = 90). The angular momentum is fixed to jz = 1/2. 
From the numerical point of view, several parameters have to be fixed. We recall that N^^r) denote the 
number of elements in each coordinates, and N* the total number of basis functions. In Fig. we illustrate 
the ground state (with i? = 1) in the prolate spheroidal coordinates, with = 4 and only N* = 20 

basis functions. Note in particular, the positions of the nuclei, at the left corners of the grid. Fig. [^reports 
the grid structure for 32 x 32 grid, and the corresponding numerical ground state. 

6.1. Convergence for TIDE 

In this section, we investigate the numerical convergence of the atomic&kinetic balance technique with 
a B-spline basis. The tests are similar to those presented in (SSj. More specifically, we study and calculate 
the ground state of Th^^®’*' (dithorium) for which Zi ^2 = 90, and HJ (dihydrogen) for which Zi ^2 = 1- The 
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CPU time / Processors: N,=N =10,20,30 



Figure 1: CPU-time / processors for time evolution of molecule (10"^ iterations) for = Nrj = 10 and = Nrj = 20 

and = Nr/ = 30 


ground state 



1 1.5 2 2.5 


I 

Figure 2: ground state, represented on 1024 X 1024 grid points 


semi inter atomic distance is set to i? = « 0.011111 a.u. for and to i? « 1.0 a.u. for H^, while 

the angular momentum is taken as jz = 1/2. The results for the calculation of the ground state binding 
energy using B-splines of order 7 and different mesh sizes are shown in Table El and for and Th^^®’*’ 
respectively. 

The results presented in this table show the convergence of the method as the number of elements, 
at fixed B-spline order, are increased. The results obtained are very accurate, although there is a small 
difference (« 10“®% and « 10“^% for and respectively) between our results and the results 

presented in [39]. This difference can be explained by a different choice of boundary conditions, different 
element formulation and different treatment of the Coulomb singularity. The B-spline basis functions, being 
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ground state: mesh 



1.2 1.4 1.6 1.8 2 2.2 2.4 


Figure 3: Mesh 32 x 32 grid points in prolate spheroidale coordinates, and corresponding ground state computed with 20 
basis functions 


polynomial with integer powers, are unable to reproduce exactly this feature. Moreover, we have that 

(79) 

(80) 


7 h « 0.999947 and V' ~ 

7Th « 0.568664 and ip ~ r-°-^3i336 


where 7H,Th are the gamma associated with a hydrogen or thorium atom. It is clear from this that the 
behavior of the wavefunction is much closer to a power law for dihydrogen and therefore, is better reproduced 
by the B-splines and thus, has a faster convergence. 

One possible cure to this is to use another prefactor in the basis function that mimics the correct behavior. 
For instance, it has been proposed to multiply the basis functions in (24) by [36l|40l[39] by 


G"(C,?7) = O 


-1+7i„-1+72 


(81) 


with 


ri = {^ + v)R, r2 = {^- ii)R. 


where 


7l:2 = 



(82) 


(83) 


and ri 2 are the internuclear distances between 1 and 2. In ground state calculations, we have = 1/2 and 
thus, 0 < 7 i _2 < 1 for Z 12 < 137. Therefore, the wavefunction has a non-integer power-law behavior close 
to the singularity at r = 0. 

The main issue with this method is that the derivative in the functionals become singular. To cope with 
this, a singular coordinate transformation can be performed that allows to transform the singular non-integer 
behavior near the nuclei to a polynomial approximation [361140] . 

From Table we can deduce the rate of convergence to the groundstate energy of reference (computed 
with iVj = Njj = 30). We represent in Fig. the logarithm of the relative error as a function iVj x in 
both cases. This graph also illustrates, that the strength of the singularity is responsible for a deterioration 
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Table 1: Results of the numerical computation for the ground state of for different mesh sizes and B-spline of order 7. Here, 
are the number of elements in each coordinates while N* is the total number of basis functions utilized. The maximum 
coordinate was fixed to ^max = 30 a.u. and the angular momentum to jz = 1/2. The calculations are to be compared with the 
results from 1391 where the authors obtained E + = -1.10264158103 a.u.. 




N* 

Min-max 

Eh+ (a.u.) 
Kinetic 

Atomic 

8 

8 

182 

-1.102590816884 

-1.102590816895 

-1.102590816899 

10 

10 

240 

-1.102638533873 

-1.102638533934 

-1.102638533914 

12 

12 

306 

-1.102641366239 

-1.102641366228 

-1.102641366222 

14 

14 

380 

-1.102641554428 

-1.102641554501 

-1.102641554498 

16 

16 

462 

-1.102641577089 

-1.102641577085 

-1.102641577079 

18 

18 

552 

-1.102641580210 

-1.102641580229 

-1.102641580219 

20 

20 

650 

-1.102641580782 

-1.102641580825 

-1.102641580823 


Table 2: Results of the numerical computation for the ground state of Th^^^"^ for different mesh sizes and B-spline of order 
7. Here, the number of elements in each coordinates while N* is the total number of basis functions utilized. The 

maximum coordinate was fixed to ^max = 15 a.u. and the angular momentum to jz = 1/2. The calculations are to be compared 
with the results from |39| and |41| where the authors obtained E + = -9504.756746922 a.u. and E + = -9504.752 a.u.. 

-*-^ 179-*-^179 




N* 

Min-max 

Etjj 179+ (a.u.) 
Kinetic 

Atomic 

8 

8 

182 

-9503.998584802 

-9504.592903093489 

-9503.999825720 

10 

10 

240 

-9504.333585765 

-9504.687718599949 

-9504.333923392 

12 

12 

306 

-9504.466070634 

-9504.711184750768 

-9504.466246166 

14 

14 

380 

-9504.539502492 

-9504.722872750701 

-9504.539637808 

16 

16 

462 

-9504.586247153 

-9504.730120406488 

-9504.586369144 

18 

18 

552 

-9504.618392312 

-9504.735095027911 

-9504.618508491 

20 

20 

650 

-9504.641636959 

-9504.738703758736 

-9504.641750168 

24 

24 

870 

-9504.672557123 

-9504.743524797539 

-9504.672667124 

30 

30 

1260 

-9504.698874401 

-9504.747650405050 

-9504.698989287 



Figure 4: Logarithm of relative error of the groundstate energy as function of x Nrj (total number of nodes) at B-spline 
order 7, for Th 2 ’’®+ and 
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of the overall convergence rate of the method for large Z. We observe in the atomic balance case, that the 
overall convergence behavior for Th^^®^ is quite similar to the min-max approach and is not as good as in the 
kinetic balance case. Although, we do not have a clear explanation for that, we think that the discrepancy 
in the convergence, is due to the presence of Vc^ with large (for small ^ 1 , 2 ; see Tablethe convergence 
rate is roughly similar to the kinetic case), in the variational intregals defining C, D, S. It is then challenging 
to numerically maintain an high order of accuracy close to the potential singularities. This will be subject 
to future investigation. 

6.2. Convergence with increasing B-spline order 

In the following test, we compute the following error on the total density, for Th^^®’*': 

pip) — II„ _ fl(p)|| 

^ ■— \\Pa P IIl2(r3_r) 

for different orders p, where denotes the numerical density constructed from an order p, B-spline function 
basis, with grid points in directions ry, and N* basis functions, and pg is a solution of reference 
constructed with very high order B-splines. In other words, we compute for fixed mesh, the error as a 
function of the B-spline order. This unusual way to show the convergence of the method is justified by 
the non-nestedness of meshes for {A(t/2®, A^^/2b j ^ 1} in prolate spheroidal coordinates, making it hard 
(without using very high order interpolation methods), to determine numerically the order of the overall 
scheme. In addition, the boundary conditions, the singularity, the number of knots, the integration method 
and its order, have all an effect on the overall order of the scheme. We here then show that the higher the 
B-spline order, the smaller the relative error, justifying the use of high order B-splines. For HJ, we report in 
Fig. ^the semilogscale of the L^-error of the overall density p{t,^,r]) = for different 

B-spIme order, and for = 6 and = Ng = 10. Results are shown using the kinetic balance 

operator. These tests show that as expected, the norm error (with respect to a a solution of reference) 
is function of the power of the order of the B-spline (semilogscale in ordinates is used). Note that in all the 
computations, the number of Gauss-Legendre points has been fixed to 64. 


Density convergence 



Figure 5: Logarithm of L^-norm error on total density as function of B-spline order (kinetic balance operator): = Ng = 10 

and = Ng = 6 


6 .3. Energy spectra of diatomic molecules 

Energy spectra are calculated using a mesh of 30x30 elements. The other parameters are set to the 
same values as in the last section where the convergence of the ground state was discussed. The value of 
the binding energies in the mass gap ([—mc^, mc^]) which corresponds to bound states are shifted by mc^ 
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to have a comparison with non-relativistic results. The values in the continua however are not shifted and 
calculated with the Rayleigh-Ritz method only. The results of the dithorium spectrum can be compared to 
the ones in [52]. Both are generally in good agreement, although a small discrepancy can be seen for the 
higher excited states. 

In the Rayleigh-Ritz method, the ribinding bound state energies shown in Tables and correspond to 
the 2N -f 1 to 2N -I- 1 -I- ribinding eigenvalues of the matrix C (once the eigenvalues are in increasing order). 
The other eigenvalues can be associated to the “discretized” negative (the first to the 27V’th eigenvalues) 
and positive (the 2N + 2-\- ribinding’th to the 4A^’th eigenvalues) energy continua. 

The convergence of the excited states is very similar to the ground state: all the values are approached 
from above and the order of convergence is close to the one of the ground state. The same is true for 
the states in the positive energy continuum, that is for E ^ mc^. For the negative energy states, the 
convergence occurs from below, but otherwise, follows the same trends as the other cases. The energy values 
in the continuua (especially their smallest and largest eigenvalues) depend on the size of the domain. In 
the dithorium calculation, the domain was smaller which yielded less accurate value in the continuua (not 
shown in the table) but better accuracy of the bound states. In all cases, the eigenvalues of the positive and 
negative energy continua accumulate at the points mc^ and —mc^, respectively. 


Table 3: Results of the numerical computation for the spectrum of for a mesh size of 30x30 and B-spline of order 7. The 
states of the positive and negative continua are computed with the Rayleigh-Ritz, Min-Max and Atomic Balance methods. The 


first 5 states are shown. 


Bound 

states 

1 

2 

3 

4 

5 


Binding energy (a.u.) 


Min-max 

-1.1026413662 

-0.6675525594 

-0.4287795568 

-0.3608697621 

-0.2554175614 


RR 

-1.1026415808 

-0.6675527718 

-0.4287811584 

-0.3608710695 

-0.2554197033 


Atomic 

-1.1026415808 

-0.6675527718 

-0.4287810919 

-0.3608690590 

-0.2553343110 


Negative 
continuum (a.u.) 

1 -18778.95240 

2 -18778.95792 

3 -18778.96471 

4 -18778.97284 

5 -18778.98233 


Positive 

continuum (a.u.) 
18778.86549 

18778.86561 

18778.86562 
18778.86741 
18778.86746 


Table 4: Results of the numerical computation for the spectrum of Th^^®"*". The mesh size is indicated on the second line. The 
B-splines are of order 7. 


States 

Naive RR 


RR 

Min-max 

Atomic 


14 X 14 

30 X 30 

30 X 30 

16 X 16 

30 X 30 

1 

-9504.6525442 

-9504.7243225 

-9504.7475523 

-9504.5862992 

-9504.6416456 

2 

-6815.3652913 

-6815.4657298 

-6815.5599111 

-6815.3230307 

-6815.3865298 

3 

-4127.8799531 

-4127.8877478 

-4128.1451137 

-4127.8197047 

-4127.8457787 

4 

-3374.4958326 

-3374.5117016 

-3374.5143753 

-3374.4569981 

-3374.4767336 

5 

-2564.1326367 

-2564.1559253 

-2564.1719708 

-2564.0744037 

-2564.0918230 

6 

-2455.9453341 

-2455.9537953 

-2455.9600280 

-2455.8837393 

-2455.9016668 

7 

-2010.6579407 

-2010.6535604 

-2010.4321103 

-2010.4241948 

-2010.4261981 

8 

-1918.5275474 

-1918.4056980 

-1915.7178408 

-1915.6761267 

-1915.6853488 

9 

-1649.5111100 

-1649.2929148 

-1643.9543595 

-1643.9320665 

-1643.9395109 

10 

-1349.5529034 

-1344.0855870 

-1313.8071916 

-1313.7606899 

-1313.7699129 

11 

-1339.1123032 

-1333.5368147 

-1303.6850950 

-1303.6580541 

-1303.6660492 

spurious 

-1218.2113620 

-1204.6990945 




12 

-1169.3956263 

-1159.1761393 

-1089.6415827 

-1089.6356220 

-1089.6370783 

13 

-1138.5709512 

-1131.0151665 

-1084.3699127 

-1084.3519981 

-1084.3522895 

14 

-1046.2053120 

-1045.4764538 

-1028.1920826 

-1028.1912423 

-1028.1920249 

15 

-1018.4013912 

-984.5252901 

-969.6816867 

-969.64172165 

-969.6482618 


Notice that although, theoretically (see Theo 3.1), spurious states can be generated using a kinetic balance 
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operator, in the numerical tests we performed with this balance, only physical states were generated, while 
using the kinetic balance. 

6 . 4 . Numerical tests for TDDE 

One important feature of the time dependent solver is the consistency with the time independent solver 
(same grid, same space discretization). A first important test is then to show that without external field, 
the density is (almost) constant in time, as expected theoretically: 

idtf) = Hofj, {t,x) e M+ X n, 4’{0,x) = tjjo{x) 

with Hq the field-free Hamiltonian and (fo eigenfunction of Hq. The formal solution is naturally: tplt, rj) = 
exp{—iHQt) ■ and we also have p{t,^,r]) = po{^,ri) = { which is the initial 

density. We numerically check that this consistency property is satisfied, discretizing the time derivative 
with a Crank-Nicolson scheme (L^—norm preserving and order 2, in time, is expected). For 8 x 
nodes, we first compute the ground state. Then using the exact same grid and spatial discretization 

(Aj = = 8, A* = 182 and B-spline order of 7), we solve idti^ = Hofj. We then report in Fig. 

{(C,??) S [0,15] X [—1,1], \p{tf,f^,r]) — Po(^)^)l}j after 10^ time iterations, with At = 10“®, that is tf — 10 . 
The solutions are represented on a 256 x 256 grid points. This result shows the strength of the time depen¬ 
dent solver with respect to its consistency with the time independent one. 


Initial density Density error: 1000 Iterations xio''^ Density error: 10000 iterations 



Figure 6: Density comparison: |p(i/, ■) — po(‘)l s,fter 10^ and 10^ iterations computed with N* = 182 basis functions and order 
7 B-spline. Representation on 256 x 256 grid points 


As a preliminary example of application, we show here the interaction of a diatomic molecule with 
a very short and intense external electric field polarized in the ^^-direction, for t ^ 0: 

A;,{t) = Aq sin^ ^ sm{uJot) 

where Nq is a positive integer, and ujq the external field freqnency. We choose in atomic units, t/ = 1, 
Aq = 100, Aq = 2 and ojq = 0.1, and the numerical data are chosen as follows: Aj = A^ = 8 and N* = 182, 
and the B-spline order is fixed at 3. We also take At = 10“^. Note that in order to precisely describe physical 
phenomena up to the zitterbewegung much smaller time step is necessary (< 10“®). We report in Fig. 
the electron driven by the field, from one center to the another, at different times (t = 0, t = 1, t = 10, t = 20). 


7. Conclusion 

This paper was devoted to the derivation and analysis of a Galerkin method using atomically or kinetically 
B-spline basis for solving the Dirac equation. We perform spectrum, as well as time dependent evolution 
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Figure 7: Initial density, poi s-nd Density variation |p(t/, •) — Po(')l time t = I, t = 10, t = 20 and with Aq = 100 


computations to illustrate some of the strengths of the method, such as its high order and the consistency 
between the time independent and dependent solvers. It was shown that using a quite reduced number of 
high order B-spline basis functions, it was possible to accuratly solve the TIDE and TDDE. This is a main 
advantage compared to finite difference methods for instance, where a very large number of points are usually 
necessary for precise computations. Atomic and kinetic balance approaches were alse compared. We recalled 
that, in term of variational collapse, the atomic balance is more relevant than kinetic one. However, due 
to additional singularities, the atomic balance was shown to be harder to accurately implemented for heavy 
ions. A future work will be dedicated to the application of the method, to intense&short laser-molecule 
interactions for pair production problems. 
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